


> 



February 1, 2008 



Abstract 



o". OPTIMAL ACOUSTIC MEASUREMENTS 

o 
o. 

^ ■ Margaret Cheney * David Isaacson ^ Matti Lassas ^ 

>■ 
O 

oo 
(N 

"^ I We consider the problem of obtaining information about an inaccessible 

,J^ ' half-space from acoustic measurements made in the accessible half-space. If 

C^ . the measurements are of limited precision, some scatterers will be undetectable 

because their scattered fields are below the precision of the measuring instru- 
ment. How can we make measurements that are optimal for detecting the 
^ ■ presence of an object? In other words, what incident fields should we apply 

^ , that will result in the biggest measurements? 

There are many ways to formulate this question, depending on the mea- 
(3 ' suring instruments. In this paper we consider a formulation involving wave- 

^^ . splitting in the accessible half-space: what downgoing wave will result in an 

f— s I upgoing wave of greatest energy? 

O ' A closely related question arises in the case when we have a guess about 

r~| ■ the configuration of the inaccessible half-space. What measurements should we 

make to determine whether our guess is accurate? In this case we compare the 
^ ■ scattered field to the field computed from the guessed configuration. Again we 

look for the incident field that results in the greatest energy difference. 

We show that the optimal incident field can be found by an iterative pro- 



^ ' cess involving time reversal "mirrors". For band-limited incident fields and 

j^ ■ compactly supported scatterers, in the generic case this iterative process con- 

verges to a single time-harmonic field. In particular, the process automatically 
"tunes" to the best frequency. This analysis provides a theoretical foundation 
for the frequency-shifting and pulse-broadening observed in certain computa- 



tions y] and time-reversal experiments |14] [15|. 
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1 Introduction 

This paper is motivated by the question "What is the best way to do acoustic imag- 
ing?" If we want to make the best possible images, we must begin with data that 
contain the most possible information. In particular, since all practical measurements 
are of limited precision, some scatterers may be undetectable because their scattered 
fields are below the precision of the measuring instrument: our data will contain no 
information about them. What incident fields are "best", in the sense that their 
scattered fields give the biggest measurements? 

This paper considers only the problem of detecting the presence of an object (or 
distinguishing it from a guess) and not the problem of making an image of that 
object. For imaging, there are other criteria for "best" that one could imagine using. 
A Bayesian criterion [^, for example, would be to look for the measurement 



producing the "narrowest" posterior distribution for the scatterer when an a priori 
distribution for the scatterer is given. 

The detection or distinguishability problem has been studied for fixed-frequency 



problems in electrical impedance tomography and acoustic scattering [|1^. The 



connection between optimal measurements and iterative time-reversal experiments 



was pointed out in |jT2[, [Tj], and [jT5[; in all these papers, the analysis was carried 
out at a single fixed frequency. The issue of optimal time-dependent waveforms in a 
special 1 + 1 - dimensional case was studied in P], where a time-harmonic waveform 
was found to be optimal. 

In this paper we study the question of optimal time-dependent waveforms in the 

3 + 1 - dimensional case. In particular, we consider the half-space geometry: we 
imagine that a plane divides space into accessible and inaccessible regions, and we 
assume that we can make measurements everywhere on the plane. 

Section 2 contains a careful formulation of the idealized problem: the wave equa- 
tion, the measurements, the notion of "biggest" . Section 3 is devoted to the example 
of a one- dimensional medium, in which the problem can be solved explicitly. Section 

4 gives an iterative experimental method that can be used to find the optimal field 
even if the scatterer is unknown. This method is precisely the iterative time-reversal 



procedure of |]T4| and [|T^. Section 5 discusses implications and open questions. The 
paper concludes with three appendices containing the technical details needed for 
the proof of convergence of the time-reversal iterates. We show that in general the 
iterates converge to a time-harmonic field that is "tuned" to the best frequency. 

2 Basic Concepts 



2.1 Distinguishability 

For any two operators Ai and A2, we say that Ai is distinguishable from A2 with 
measurement precision e if the distinguishabihty 6{Ai,A2), defined as 

d{Ai,A2) = snp :-— , (1, 

/ 



is greater than e. A field that is best for distinguishing Ai from A2 is an / for which 
the maximum is attained. We will determine below the norms that are appropriate 
to use. 

2.2 Acoustic Wave Equation 

We consider the constant-density acoustic wave equation 

{V^~c-\x)d^)U{t,x) = 0. (2) 

in the case in which c = Cq everywhere in the upper half-space X3 > 0. This model 
includes neither dispersion nor dissipation. 

We can formulate the scattering problem in a variety of ways 0. In particular, 
we can use either a boundary map, sources, or a scattering operator defined in terms 
of wave splitting. 

2.3 The Boundary Map 

To define the boundary map, we specify that U = f on the surface X3 = 0. This condi- 
tion, together with an outgoing radiation condition at infinity P], uniquely determines 
a solution U in the lower half-space. We can then take the normal derivative dU/dx^; 
this normal derivative, restricted to the surface X3 = 0, we denote by g. The mapping 
from f to g is the boundary map A. Thus, on the surface X3 = 0, AU = dU/dx^. 
Note that A is an operator-valued function of time. 

Acoustic distinguishability can be defined in terms of the boundary map as 

||(A-Ao)/|| 
dB{c, c") = sup f^-p^^ (3) 

/ 



for appropriate norms. Here Aq denotes the boundary map for the reference sound 
speed c°(x). This formulation, in terms of the boundary map, is not pursued in this 
paper. 



2.4 Sources 

To formulate scattering in terms of sources, we consider the wave equation with a 
source: 

{V^ - c-\x)dl)Uj{t,x) = J{t,x). (4) 

Scattering data is then Uj{t, x) for x on the plane, where J is supported on or above 
the plane. Acoustic distinguishability in terms of sources would be 

(55(c,c°) = sup l'^;~^°" , (5) 

J IkII 

where C/j denotes the field due to the reference sound speed c^{x) and source J. This 
formulation is not pursued in this paper; instead we consider the scattering operator. 

2.5 The Scattering Operator 

We define the scattering operator in terms of upgoing and downgoing waves. The 
motivation for this point of view is the existence of network analyzers, which can de- 
compose a time-harmonic signal in a waveguide into an upgoing one and a downgoing 
one, and measure the amplitude and phase of the upgoing wave. Stepped-frequency 
radar, for example, is based on the ability of such instruments to transmit and receive 
signals at the same time. 

2.5.1 Upgoing and downgoing waves 

To define upgoing and downgoing waves, we make use of two Fourier transforms, a 
temporal one and a spatial one. First we inverse-Fourier transform the solution U of 
d) in t: 

u{u, x) = J^-^U = {2tt)-^ I U{t, x)e''^'dt. (6) 

This frequency- domain solution u satisfies the reduced wave equation 

(V^ + cjV2)m(cj,x) = 0. (7) 

We write k = uj/cq and with a small abuse of notation we write u{k, x) instead of 
u{cok,x). We then Fourier transform u again in x' = (2:1,0:2), so that 

u{k, 7]', X3) = F^^u = f u{k, x)e-^'='''-^'rfV (8) 



where t]' = {r]i,ri2). (We note that F^.^ depends also on k.) Then U is recovered as 
U{t, x) = -^ J J u{k, T]', X3)e"''^'-^'e-'''^'ed^rj'codk. (9) 



In the upper half-space, it satisfies the ordinary differential equation 

(dl + k'-k'\vf)u = 0, (10) 

which has the general solution 

u{k, 7]', xs) = A{k, riy^'^^''^ + B{k, V)e-^^''^^^ (11) 

where 



, ,, 1 — |r7'P for 1 > \rl\ , ^ 

^(sgn/cjw |?7'p — 1 tor 1 < \ri \ 



We define the vectors r]^ = (r/', ±773), which satisfy r]^ ■ t]^ = 1. 

In order for U, as defined by (^, to be real- valued, the Fourier transform u must 
satisfy certain symmetry conditions. In particular, we must have A{—k, rj') = A{k, t]'), 
B{-k,r]') = B{k,T]'), and r]s{-k) = rjs^k). 

Equation (pUf ) shows us how to split the time-domain solution of (|]) into two 
parts, which we call the upgoing and downgoing parts. Thus for 2:3 > we write 
U = W + U^, where U^ is 



U\t,x) = A{k,r]')e"''^ ■''e-''"'"'k^d^7]'codk, (13) 

and 

f/i(t, x)= I I B{k, r/Oe'^^'-^e-^^^o^PrfVcoC^fc- (14) 

We see that (jl^) and (jl^) are plane wave decompositions. The components for 
which \ri'\ < 1 are propagating plane waves, and rj^ is a unit vector that gives the 
direction of propagation. The sign of the third component of r]^ determines whether 
the wave is downgoing or upgoing. On the other hand, components with \ri'\ > 1 cor- 
respond to evanescent waves. For U^, these evanescent waves decay in the downward 
(negative 0:3) direction; for U\ they decay in the upward (positive) direction. 

2.5.2 The scattering operator 

It is natural to define a scattering operator S^ as the map from U^ to U^ . We denote 
the kernel of this operator also by S^: 

U\x,t) = r fsHx,y,t-T)U\y,T)d^ydT. (15) 

J —00 J 

We note that this scattering operator is defined only on downgoing solutions of the 
Helmholtz equation, i.e., on functions of the form (|1^). 



The kernel of (|I5|) is a convolution in time because the Fourier transform "diag- 
onalizes" the time derivative of @, so that the frequency is simply a parameter in 
(|^). In other words, the convolution is an expression of the fact that in the frequency 
domain, ([T5|) takes the form 

u^{k^x)= f S^{k,x,y)u\k,y)d?y. (16) 



The time-domain operator S^ is related to the frequency-domain scattering operator 
Si by Si = T-^Slj^. 

In fact, Sl and Sl are determined by their actions on the plane X3 = 0. We see this 
as follows. First we Fourier transform ( [I6| ) in space. The operator Sl is transformed 
into the operator S^^^stz = Fx^Sl-F^.^ = {Fx^J^~^)Sl (F^^^J^^^)^^ , where F^^ is defined 
by dj) and JF~^ by (^). The transformed version of (^) is 

u\k,r]',X3) = I S^^^s:.^{k,r]',fj')u^{k,fj',X3)k^d'^fi'. (17) 

The operator Sx3,x3 is determined completely by its action at X3 = 0, which we see 
from the following argument. 

Into ([T7|) we substitute «-'■(/;;, 77', £3) = B{k,fi')e~'''^^^^^ and u''{k,ri',x-i) = A{k,r]')e^^^^^'^] 
we see that 

S{k,i,fi') = e-^'=''^^^S',3,,3(A;, V,77')e-'^^^*« (18) 

satisfies 

A{k, 7]') = [ Sik, ri\ fi')B{k, fi')Pd^fi' (19) 



where A and B are as in (|ll|). The relation between the operators S and Sx3,x-i can be 
written S^-^^x^ = Ex^SE^^, where E^^ is the operator of multiplication by exp(2/i;?73X3), 
and thus 

Sl = F'^Ex^SEx^Fx^. (20) 

On the plane 0:3 = 0, this becomes 

S = F,-'SFo (21) 

This defines a scattering operator S on the plane x^ = 0. It is this operator, together 
with the corresponding time-domain operator S = J-'SJ-'~^, that we will use in the 
rest of the paper. We note that the domain of the operator S is restricted to the 
space of downgoing waves as defined by ([I^) . 

The scattering operator S and the boundary map A are related to each other by 
formulas developed in (See Appendix A for details). They are thus equivalent 
operators. Whether it is better to formulate a given problem in terms of a scattering 
operator or a boundary map depends largely on the design of the equipment involved. 
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2.6 The Energy Identity and the Energy Flux 

If we multiply (H) by dtU and integrate the resulting equation over the volume V, we 
obtain 



We write the first term of (P) as V ■ {{dtU)VU) - V{dtU) ■ Vt/, and apply the 



divergence theorem to the term containing the divergence. We thus obtain 

l/,_..,o 1 



/ {dtU)d. 

JdV 



WdS = dtj^-{\VU\' + —-(dtUfyx, (23) 



where u denotes the outward unit normal to the surface dV . 

This equation relates the change in energy in the volume V (the right side of (p3D ) 
to the energy flux across its boundary surface dV . 

From (^3]) we see that the time-integrated energy flux across a surface dV in the 
normal direction v is 

/oo /• 
/ {dtU)d^UdSdt. (24) 

-ooJdV 

We can use Parseval's identity to write the time-integrated energy flux in terms of 
the frequency-domain wave functions: 

/ iiuu)d^u dSduo = -(2nf / / (icoku)duU dScodk, (25) 

-oo JdV J-oo JdV 

where the overline denotes the complex conjugate. 

2.6.1 The energy flux of upgoing and downgoing waves 

At the surface X3 = 0, the total field splits into upgoing and downgoing parts. 
Because the flux is quadratic, it does not necessarily split into corresponding up- 
going and downgoing fluxes. However, a quick calculation using (|T^, ([T^)? and 



([25| ) shows that if upgoing and downgoing evanescent waves are not both present 
on the plane X3 = 0, the time-integrated cross terms j'^^ j^^^Q{dtU^)dx.JJ^ dx' dt and 
/-^ Ix3=oi^tU^)^x3U^dx'dt cancel. Under these conditions, the time-integrated fluxes 
do split into upgoing and downgoing fluxes, so that W{U^ + W) = W{U^) + W{W). 
Throughout this paper we assume that the sources of the downgoing fleld are far 
from the scatterers, so that there is no interaction between upgoing and downgoing 
evanescent waves on the plane X3 = 0. The flux of a downgoing wave is positive; that 
of an upgoing wave is negative. 



We write the downgoing energy flux as 



oo 



W{U^)= / {dtU^)d^,U^dx'dt ={27rf \B{k,r]')\^clr]3k^d^r]'dk, (26) 

where we have used ( pjlf ) and (^51) in carrying out the computation (^61) . In (|26|) there 
is no minus sign because downgoing energy travels in the — X3 direction. 

Although the left side of (p6|) is real, it is not obvious that the right side is, because 
?73 can be imaginary. However, if one splits the k integral into pieces as 

W{U^)=m' + j^^\B{k,i)\''clk%dkd^i (27) 

and uses the symmetry properties of B and 773, one sees that for evanescent waves, 
the two terms cancel. This shows that the evanescent waves do not contribute to the 
energy flux. 

The flux can be used to form an inner product on the space of downgoing propa- 
gating waves; we define 

I J -00 Jx3=0 ^ ' 

= (27r)3 / / u^¥clk%dkdS (28) 

J J\ri'\<'^ 

We note that the product clk'^rj^ is non-negative, so in the transform domain, this 
inner product is merely a weighted L^ inner product. 

Similarly, the energy flux of the upgoing scattered field SU^ = W that passes 
through the plane X3 = is 

/oo r yoo r 

/ idtW)d^^Wdx'dt = -{27rf \A{k,ri')\^clk%d^ri'dk. (29) 

-00 Jx3=0 J —OD J 

The minus sign in (^9]) is due to the fact that the upgoing wave corresponds to energy 
leaving the lower half-space. Equations (p6D and (p9D show that the time-integrated 
downgoing flux W{U^) is positive and the time-integrated upgoing flux W{W) is 
negative. 

Note that for propagating waves, \W\ satisifes the triangle inequality: |W^(6^i + 

f/J)|<|W^(f/i^)| + |W^(f/J)|. 

The flux inner product on the space of upgoing propagating waves is 



d X dt 



I J-00 JX3=0 ^ ' 

= (27r)3 / / uWclk%dkd^7]'; (30) 

thus for propagating waves, |Vr(f/^)| = {W,W)w 



Conservation of energy. If the medium is initially quiescent, conservation of en- 
ergy tells us that W{U^) > \W{SU^)\; this can be seen from integrating ( ^3]) over all 
time, and using the fact that the energy within the volume V is initially zero, and 
cannot become negative. The time integral of the right side of (^) is thus positive. 
The left side we write as W{U) = W{U^) + W{SU^). This implies that the total 
upgoing flux |iy(iSf/^)| cannot be greater than the total downgoing flux W{U^). 

Finite-energy fields on the plane. We deflne the space w of finite-energy func- 
tions on the plane to be the closure of C^(R^ x R) in the inner product 



(u, v)l = j j u{k, r])v{k, r])clk'^\r]3\d'^T]'dk, (31) 

and the space W = J-'w. 

2.7 Acoustic Distinguishability via the Scattering Operator 

We define the acoustic distinguishability in terms of the upgoing and downgoing 
energy fiuxes through the surface X3 = 0. 

For a reference scatterer with scattering operator Sq, the energy fiux of the upgoing 
scattered field SqU^ and of the difference field {S — Sq)U^ are defined similarly. 

In general the distinguishability of S from Sq with the incident field U^ is 

,(5.5.)^supM(£_^.,„pM(|_^. ,32) 

We recall that evanescent components do not contribute to the energy fiux. To 
remove the evanescent components from (^2]), we denote by P the orthogonal pro- 
jection onto the propagating components: P = Fq^PFq, where P is the operator 
of multiplication by X\ri'\<i, the function that is one for \t]'\ < 1 and zero otherwise. 
Explicitly, P is given by 

P{k)f{x') = fk' f e''^'<^^'-y'^d\f{y')d'y'. (33) 

In the time domain, V = J-'PJ-'^^. With this notation, we can write W{U) = 
W{rU) = W{Pu). 

Moreover, the scattered field due to an evanescent incident wave has zero total 
energy fiux. This is because of the comments at the end of the previous section: 
= W{{I - P)u^) > \W{S{I - P)u^)\ implies that W{S{I - P)u^) = 0. 

In addition, the upgoing energy fiux can only be increased by getting rid of the 
evanescent components of the incident wave. This is because of the triangle inequality 
\W{SPu + S{I - P)u)\ < \W{SPu)\ + \W{S{I - P)u\ = \W{SPu)\. 
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This implies that the downgoing waves that give rise to the maximum total energy 
flux are propagating waves. Thus we find that the distinguishability can be written 

\W{V{S~S,)VU^)\ _ \W{P{S-S,)Pu^)\ 
5{S, So) - sup ^^^^;^ - sup ^^^^ . (34) 

We note that the scattering operator Sq for free space is the zero operator. Thus, 
according to (0) and (|32|) , the presence of a scatterer can be detected with measure- 
ment precision e if the distinguishability satisfies 

\W{VSVU^)\ \W{VW)\ ^ 

^^^' °) = T wivui) = J",p mvmy > " ^''^ 

The distinguishability can be defined equally well in terms of the operator S of 

morsoim- 



3 Example: The One-Dimensional Case 

If the medium in the lower half-space depends only on depth, then the coefficient B 
of ( ]TT| ) is the refiection coefficient R{k, t]') multiplied by the incident coefficient A. In 
this case, the distinguishability (5(5,0) can be computed from (p6D, (p9D, and (|35| ) as 

^^'^'"^"T JT\Bik;wWmdkd^' ■ ^ ^ 

The maximum of the right side of (^) is attained in the limit when i? is a delta 
function supported at the maximum of \R\. 

Thus to maximize the scattering from a one-dimensional scatterer, we compute 
the conventional refiection coefficient R{k, rj'), and find the values of k and r]' at which 
it attains its maximum. Taking i? to be a delta function supported at these points 
corresponds to taking an incident field that is a plane wave of fixed frequency uj = Cok 
and incident direction given by rj'. 

Note that since R is minus one for \ri'\ = 1 (grazing), a maximum always occurs 
at grazing incidence. If this is undesirable, grazing incidence can be excluded by 
modifying the definition of distinguishability. 



4 An Adaptive Method for Producing the Best 
Fields 

To maximize the distinguishability when the medium is unknown, we can use the 
following adaptive method. 
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We write 
jr,c:c,_ \W{V{S-So)VU^)\ _ {U, {V{S - S,)VyV{S - So)VU)w 

^^^' ^'^ - T mvm) - S (vuTpUw 

(37) 
where U{x,t) = U^\x.^=o and where the adjoint * has been taken in the space W 
(defined just below (pl|)). 

We see in Appendix A that (VSV)* = T(VSV)T, where T denotes the time- 
reversal operator TU(t,x) = U{—t,x). 

Thus we see that the operator appearing on the right side of (|37|) is ^ = TV{S — 
Sq)VT{S — So). In general, to maximize a quotient of the form 

{U,AU)/{U,U), (38) 

one considers an appropriately normalized sequence A"'U. When A is compact, this 
sequence converges to the largest eigenvalue of A. Here, however, A has a contin- 
uous spectrum, so we expect the sequence A^U, when appropriately normalized, to 
converge to a generalized eigenfunction of A, and the corresponding quotient ( p8| ) to 
converge to the supremum of the continuous spectrum. We note that such generalized 
eigenfunctions do not have finite energy. 

When A is compact, the usual way to normalize A^U is to divide by H^^f/H . Here, 
however, because we expect A^U to converge to a distribution in the time variable, 
we must use a distributional normalization. We consider test functions in a particular 
space that is discussed in the appendix. These test functions are functions of space 
and time. In the time variable, they are Fourier transforms of functions of compact 
support. The distribution action is chosen to coincide with the fiux inner product 
defined by (pSf ) and (^). For the distribution action we use the same notation {■,-)w 
as for the fiux inner product. 

To normalize, we choose an arbitrary test function \E', and consider the sequence 
A"'U/{A"'U,'^)w This gives rise to the following algorithm for carrying out the 
maximization of (pT]). 

1. Start with any Vq] let j = 0. 

2. Send Vj into the lower half-space, and measure the resulting upgoing field 
Vj{t,x) =SVl{t,x). 

3. Calculate the corresponding scattering from the reference configuration SqVj {t, x)). 

Calculate the difference field Vj{t,x) = Vj{t,x) — SoVj'{t,x). 

4. If j is even, let 

Vl+,{t,x) = Vj{-t,x), (39) 

add one to j, and go to Step 2. 
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5. If j is odd, normalize: 

''-<*-' -US' '^°' 

add one to j, and go to Step 2. 

Appendix B contains a proof that, in the case of a compactly supported scatterer 
in free space, the sequence f/„ = V2n = A^U generally converges to a single time- 
harmonic wave. The frequency of this wave is the frequency at which the largest 
eigenvalue of 5* attains its maximum. If this largest eigenvalue happens to attain 
the same maximum at several different frequencies, then the iterates Un converge to 
a sum of time-harmonic waves with these frequencies. The relative strengths of the 
different frequencies is determined by the corresponding frequency components of the 
initial incident wave Uq = Vq. 

The argument in Appendix B takes place within a limited frequency band; this 
frequency band is determined by the bandwidth of the test function. 

We note that as expected, the limiting time-harmonic waves do not have finite 
energy. This property also appears in the one- dimensional example (|36|) . 

Step 4 can be omited and Step 5 performed for every j: the linearity of the problem 
implies that extra normalizations do not affect the limit. The proof in the Appendix, 
however, corresponds to the above algorithm. 

The algorithm can also be implemented including a step in which the evanescent 
waves are filtered out. If they are not filtered out, however, they will die out anyway 
as the iteration proceeds, because experimental time-reversal of a field that includes 
evanescent waves is simply another physical field with evanescent waves. 

5 Conclusions and Open Questions 



This analysis shows that the iterative time-reversal work of ||T4| and [|I^ provides an 
experimental method to obtain optimal fields. Moreover, this analysis explains the 
frequency-shifting and pulse-broadening seen in |]TB[ and [^ : the optimal time-domain 
waveform is a time-harmonic one tuned to the best frequency. 

This analysis suggests that the commonly-used pings and chirps are not optimal 
from the point of view of distinguishability. 

There are many open questions related to this work, one of which is the question 
of limited-aperture and limited-time measurements. Upgoing and downgoing waves 
in a limited aperture can be defined with the help of eigenvalues of the Laplacian for 
the aperture. However, it is not clear how to determine the entire incident wave if the 
incident wave is known in only a limited aperture. This involves a detailed modeling 
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of the transducer or antenna. Perhaps a formulation in terms of sources will be more 
useful in this case. 

We have not studied the question of whether the distinguishability, as a function 
of the medium, is monotone in any sense. This is an important issue for the following 
reason. Suppose we discover that a sphere of a certain radius is detectable with 
a certain measurement precision. Does this imply that a larger object will also be 
detectable? For fixed-frequency measurements, the answer to this question is certainly 
no, because of the phenomenon of resonance. A small sphere may happen to have 
a radius commensurate with the wavelength of the probing wave, and may therefore 
scatter much more strongly than a larger sphere. However, the use of time-dependent 
fields may give different results. 

The simple wave equation studied in this paper does not include the important 
effects of variable density, dispersion, and dissipation. 

Moreover, the question of distinguishability is only the first step in building an 
optimal imaging system. How should we choose a full set of optimal fields that could 
be used to form an image? 
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A Appendix: Properties of S. 

A.l Expression for Kernel 

We can find an expression for the kernel S of ([T7|) and (p!8[) by taking u^ to be a 
delta function. The kernel of S is then the corresponding upgoing wave u\ Taking 



u 



to be a delta function means that we take u^ of the form exp(iA;?7~ ■ x) for some 
V~ = {v'-i~V^)- Here fj^ can be complex. We write the corresponding frequency- 
domain field as ip: 

iPik, X, fi-) = — j U{t, X, fi-)e''""''dt (41) 

so that 

u{k, V, X3, r/-) = / ^{k, X, fi-)e~'^'''-^'£x'. (42) 
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In the case of scattering from a perturbation in free space, we can express the 
total frequency-domain field ip as a solution of the Lippmann-Schwinger equation 



ip{k, X, fi) = exp{ikr] ■ x) - k J g{k, x, y)V{y)i){k, y, r])d y, (43) 

where g is the usual outgoing Green's function 

^ik\x-y\ 

g{k,x,y) = — (44) 

47r|x — y| 

and V{y) = 1 — cl/c^{y). The scattered field u^ thus is represented by the integral 
term of (||). 

The Green's function can be written in terms of its two-dimensional Fourier trans- 
form as 

g{k, X, y) = — ^ / _J_^^krJ3\x,-y,\^ikrJ'■{x'-y')J^2^2^>_ ^45^ 



To compute A and B of (pH]), we take the xi, X2 Fourier transform of (|43| ) in the 
region X3 > 0. We assume that the perturbation V is supported in the region ^3 < 0, 
so that when we consider (^) we can remove the absolute values in (^). In the 
transform domain, the scattered field is given by 

u\k,r,\x^,r) = -k' j ^J e-''^^-yV{y),lj{Ky,r)dydy^e''^^''\ (46) 

This shows that 

5(fc,V,r/') = -T^A(A:,r^+,r/-) (47) 

2r/3 

where 

A(fc, 77, 77) = I e-'''^-yV{y)ij{k, y, ^)Sy (48) 

is a scalar multiple of the classical scattering amplitude ||13|. (This A is not to be 
confused with the A of (|TT])!) It satisfies the reciprocity relation 

A{k,ri,f])=A(k,-f],-ri) (49) 

and the symmetry relation (for real-valued perturbations V^ 



A{k,i],f])=A{-k,"q,f]). (50) 



It is clear from (^71) that S is an analytic function of k [|1^ . 
An expression similar to (^) can be obtained for the field scattered from a per- 
turbed half-space or layered medium; in this case the appropriate background Green's 



function should be used instead of the free-space Green's function in (^31) . Scattering 
theory in such cases is considered, for example, in |T9|, [Q, and 0. 
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A. 2 The Adjoint 

For a scatterer in free space, the adjoint of VSV in the space W can be computed 



exphcitly as follows. From (|30|) we have 



oo 



{VSVU,V)w = {2ny / / S{k,r]',fi')u{k,fi')k'd'f]'v{k,7]')dk%dkd'r]' 

-ooJ\r)'\<l J\fi'\<l 



oo 



(27r)M / u{k,fj') S{k,r]',fi')v{k,7]')7]3d'^T]'clk'^dkd^fi' 

J-ooJ\rt'\<l "'|»?'|<1 

(51) 



From §3), (H), and (|0D, we have 
{VSVU,V)w 



(27r)3 J-ooJ\-n'\<i ' "'|^'|<i 



u{k, fj') / ~{:ik/2)A{k, -f}~ , -r]+)v{k, ri')d'^r]'clk'^dkd^f]' 



u{k, fj') / {ik/2)A{-k, -fj-, -r]+)v{k, r]')d^ri'clk^dkd^' 

-ooJ\ri'\<l A'n'\<'^ 

oo r 7 ~ 

/ u{k,fj') f S{—k,—fj',—7]')v{k,7]')d'^r]'rj3{—k)c^k^dkd^fj', 

'OO J|7y'|<l "'h'|<l 

(52) 

where in the last equality we have used the fact that S{—k, r]', fj') = ik{2ri3)^^A{—k, r]^, f)^). 
We note that in the course of this computation, the r^s has disappeared and has been 
replaced by fj^,] this is because of the 773 in the denominator of (^7|). 



We see that the action of the adjoint {VSV)* on V is given in the transform 
domain by the rj' integral of (|52|) : 

FT-\VSV)*J^F-^v{k,fi')= f S{-k,-fi',-7]')v{k,7]')ed\. (53) 

J\v'\<i 

However, S{k, t]', fj') is defined by 



00 /-oo 



6 (co{k - k)) S{k, T]', fi')= / / / e-''='''-"'e*'^'=°*5(t-r, x', y')e-'^''^^e'^^'-y' d^'d'^y'drdt. 

(54) 
From this, we see that 



5 [coCk - k)) S{-k, -T]', -fi') = / / / e-'^^'"'e-*'=^°*5(t-r, x' , y')e'~^'°^e"''^'y'd^x'd%'dTdt. 

(55) 
Letting t —>■ —t and r —>■ —r in (155|) shows that the kernel of (|53|) corresponds to the 
operator TVSVT, given by 

/oo /"OO 
/ 5(r-t,x',y')^(^,2/')rf^A'- (56) 

-00 J —00 

Thus (P5P)* = T{VSV)T. 
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A. 3 Compactness 

Theorem 1 Assume that the sound speed c{x) is bounded and differs from cq only in 
a bounded subset of the half-space X3 < —h < 0. Then the fixed-frequency scattering 
operator S is compact on the weighted space 



Ll = {f:f{v')\l-W\r'eL'}. 



(57) 



Proof. We use (^71) and (|T^) to compute the square of the Hilbert-Schmidt norm 



of S in the space Lg 



2. 



WSfn.. 



4 



^'^■yV{y)^{k,y,fj')d' 



1-1^? 



i'\2 



1/2 



dWd'f]'. (58) 



We use the fact that ip can be spht into an incident and scattered field via 
This allows us to split the kernel ( |¥7D into two parts: 

S{k, T], fj') = Ssik, T], fj') + Ssc{k, ri, fj') 

where the "Born" term is 

ik 



and 



Ssc{k,'n,v' 



SB{k,7],T]' 



ik^ 



2r/3 



Mn-vYvy^y^^d^y 



2m 



''"^■'Viy) / giy - z)Viz)ij{k, z, fj')d'zd'y. 



(59) 



(60) 



(61) 



We compute the Hilbert-Schmidt norm of each part. 
The Hilbert-Schmidt norm of the Born term is 



ll'S's||ii-.5. = -r 



Mv-Vyyv(y)d^y 



1 — \t]' 



r^'\2 



1 — Irj- 



'12 



1/2 



d^T^'d'^fj'. 



(62) 



The y integral of ([6^ ) can be written 



^iKn' -n')-y' ^-ik{m+m)yi 



V{y)d^y. 



(63) 



For \ri'\ < 1, for which 773 is real, this integral is bounded when V is integrable. For 
\rj'\ > 1, the integral decays exponentially because V is supported in the region where 
ys "^ h < 0. The same comments apply to the behavior in fj'. Thus this integral is 
bounded by cexp{—hk{\ri'\ + \fi'\)). This estimate can easily be used to show that 
II'S'bIIh.s'. is finite. 
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The Hilbert-Schmidt norm of the scattered part is 



19 IP 

\'-'sc\\h.S. 



4 



'''■'V{y) / giy - z)V{z)^{k,z,fj')dhdh 



1-1^? 



1 — |?7' 



/|2 



1/2 



cfrj'cfri'. 
(64) 



An apphcation of the Cauchy-Schwarz inequahty shows that the y integral of 
is bounded by 



'"'■'ViymLHy) \Viy)\'/' g{y - z)V{z)^{k,z,r,')d'z 



L^y) 



(65) 



Direct computation shows that the first norm in ( |65|) is bounded by cexp(— /iA;|?7'|). 
Standard scattering theory arguments (See Appendix C) can be used to show that 
the second norm appearing in (|65|) is bounded by cexp(— /iA;|f/'|); thus the y integral 
satisfies the same bounds as in the Born term. QED 

Remark. In [||, it was shown that on L^, the operator PSP has norm less than 
or equal to one. 

For future reference, we note that the operator a{k) = (PSP)* [PSP) is given 
explicitly as 



{Fam' 



l^'j'l<i-'IC'l<i 



Si-k, -fi, -C)S{k, c, v')dXfiv')k'd'v'- 



(66) 



B Appendix: Convergence of the Iterative Algo- 
rithm 

For simplicity of notation we consider only the case when Sq = 0. In this case, the 
nth iterate is 



Ifl 



{T{VSV)T{VSV)YU^ 



A'^U^ 



i 



{{T{VSV)T{VSV)YUi ^b)w (^"f/o% "^b] 



(67) 



w 



where A = (VSV)*(VSV), where the star denotes the adjoint with respect to 
the flux inner product (•, ■)vk- We simplify the notation by dropping the arrow 
on U. Because we expect the limit to be a distribution, we consider the quantity 
{Un, $)iy = J{un, 4>)x'clk^dk, where $ is a smooth test function, = JF$, and ( , )x' 
denotes the weighted inner product that for smooth functions is (/, 0)^/ = (/, 0)^2 = 

Specifically, we consider test function that are functions oft and x'. When inverse 
Fourier transformed in t and Fourier transformed in space, at each frequency they 
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must be in L^, and in the frequency variable they must be integrable and (uniformly) 
supported in the compact interval [—B, B]. We denote this space of test functions by 
X. 

In the frequency domain, the nth iteration is J-'^^{A'^U) = a"'u, where JF~^ denotes 
the inverse Fourier transform (^) and a = T^^AT = {PSP)*{PSP). From our 
normalization of f/„, we have 

_ {A^U,<I>b)w _ J{a^u,<f)B).'clk^dk 

We note that A and a{k) are self-adjoint on the space W and on L^, respectively. 
Moreover, Theorem 1 shows that a{k) is compact on L^, and it can therefore be 
written a{k) = J2i M{k)Pi{k), where A; > A/+i and the Pi are orthogonal projections. 
Because A is non-negative, all the eigenvalues A; are non-negative. We see from (|5B|) 
that a is analytic in k, and the As are therefore piecewise analytic p. 

Suppose that Xo{k) attains its maximum in the set {\k\ < B} at ko, and that 
A(A;) = M. Then in a neighborhood of ko, X{k) has a Taylor expansion whose first 
two terms are M — b{k — k^Y for some h and some integer p. We call p the order of 

Ao. 

We allow eigenvalues with different indices to coincide at a point; thus it is possible 
that a finite number of eigenvalues also attain the maximum M at fco- Iii this case, 
these eigenvalues have a Taylor expansion similar to that of Aq, possibly with different 
fe and ps. In this case we also refer to the relevant p as the order of the eigenfunction. 

We will need the following lemma. 

Lemma 1 Assume that b is positive and that p is an integer. Then for large n, 

I{n,p) = j\l-bkTdk^^^, (69) 

where C{p) is a nonzero constant independent of n. Thus the convergence to zero of 
I{n,p) is slower for larger p and smaller b. 

Proof. Let s = b^^'^k. Then / = b^^^^ /q (1 — s^)"ds. Replacing the upper limit 
by 1 results in an error that is exponentially small in n. Denote by /„ the integral 
/q^(1 — s^)"ds. Then we can write 

In+l = f {1-SP){1- sP^ds = In- f SP{1 - S^Yds. (70) 

Jo Jo 

In the integral of (^), we integrate by parts, differentiating s and integrating s^~^(l — 
s^Y- The boundary term vanishes, and ( [TDD becomes 

In+l = In-\ 7 —rrln+l- (71) 

p{n + 1) 
18 



Solving for /„+i gives the recursion 



Since Jq = 1, we have 



^'^^^ - p(I + l) + /- ^^^^ 



^"-(pll) (2pIl)'"(npTl)- ^^^^ 



Taking reciprocals and logs and expanding, we find that 

n 
-login = ^log(l + l/(ip)) 



-E--E^Er^ (74) 



Exponentiating and taking reciprocals again, we have 



1 " 1 
J„ = C(n,p)exp(— ^-), (75) 

Pj=iJ 



where C{n,p) has the large-n limit 



/ oo (_l)m \ 

lim C(n,p) = exp - ^ MrC(^) , (76) 



m=2 "^P" 



where Cl'^) denotes the Riemann-zeta function ({m) = J2T=i k~"^. Thus we see that 
the large-n behavior of /„ is determined by the second factor of ([75|). 

We determine the large-n behavior of this second factor as follows. From approx- 
imating the sum by a Riemann integral, we have the estimate 

n ^ 

logn < log(?2 + l) < y - < 1 + logn (77) 

We multiply by —1/p and exponentiate to obtain 

1 " 1 
g-i/p^-i/p < exp(-- y -) < n-i/P (78) 

Pj=iJ 



QED 
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Theorem 2 Assume that a is an analytic self-adjoint-compact-operator-valued func- 
tion of k having the representation a{k) = J2i-^i{k)Pi{k), where Xi > A^+i and the Pi 
are orthogonal projections. Assume that Aq is not a constant function ofk. Then for 
test functions (pB^k, x') and ipsik, x') in X whose support in the frequency domain is 
in the set {\k\ < B}, 

^.^ J{a''u,(j)B)x'kdk ^ Ei,j(3ij{Piu,(j)B)x'{kj) 

n^oc /(a"M, ipB)x'kdk J2ij PljiPlU, i'B)x'ikj) 

where the sums are over those indices j and I for which Xi{kj) = M, where M is the 
maximum of Xq in the set {k : \k\ < B}, and for which Xi has maximal order at kj. 

Proof 

The representation for a allows us to write (^) as 

(TT ch \ IElXf{PlU,(pB)x'kdk 

lEiKiPi'^^i'BJx'kdk 

The A; and Pi are peicewise analytic functions of /c [^. In particular Ao(fc) is piecewise 
analytic, and therefore attains its maximum M on a discrete subset of the set {\k\ < 
B}. We cover the support of (p with open intervals Nj so that each Nj contains only 
one kj. We decompose the test function 0^ as 4>b = J2j(pj 0) where the 0j are in 
C^{Nj) and (f)j = 0^ in a neighborhood of kj. We carry out a similar decomposition 
for ipB- 

With the notation fijik) = k{Piu,(f)j)x' and gijik) = k{Piu,iljj)x', we can write 
as 

Ei,jlN,Kik)fij{k)dk 

Ei,jlNjKW9i,jWdk 

We divide the numerator and denominator of (^Tj) by M", and write ri{k) = 
Xi{k)/M] thus |ro| < 1, and |r/| < 1 for all but a finite number of values of /. Then 
([STP can be written 

Ei,jlN,rt{k)fij{k)dk 

T.i,jSNjf{k)9i,j{k)dk 

We write 

I?, = I rnk)fi,{k)dk. (83) 

J Nj 

We multiply and divide /^ by Jj^. rQ{k)dk, and write (J^ = r^ I f^ r^dk. 
For those / = 0, 1, . . . Lj for which A; attains the maximum M and thus ri attains 
the value 1, we add and subtract fi,j{kj) to the quotient, obtaining 

I?, = (^fiAkj) + I QmkAk) - fiAkj))dk^ I r^{k)dk (84) 
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We will show that the integral term within the parentheses on the right side of (p^) 
vanishes as n goes to infinity. For this we use the two facts that 1) except aX k = kj, 
Q converges to zero pointwise as n goes to infinity; and 2) J (^ = 1 for all n. 

Given e > 0, we choose N^ so small that on Nj, \fij{k) — fi,j{kj)\ < e/2. The 
integral in parentheses on the right side of (^^ we split into two integrals, namely 
A^ and -B", where 



^n 



N' 



C?{k){fi,,{k)-kAk,))dk, 



and 



B' 



NAN' 



QikKfi.ik) - f,,{k,))dk. 



15) 



(86) 



Then A^ < e/2. Next, we choose N so large that for n greater than A^, B"^ < e/2. 
This shows that the integral in parentheses on the right side of (Q) vanishes as n 
goes to infinity. Thus (|^) is a product of a factor converging to fij{kj) and a factor 
converging to zero. 



For I = Li + 1, 



Lj + 2, . . ., for which A; is strictly less than M, we write 



'i,j 



Q{k)fi,{k)dk] / rl{k)dk. 



N, 



^7) 



In this case, (^ converges to zero pointwise for all k. Thus the integral in parentheses 
of /"• converges to zero by the Lebesgue Dominated Convergence Theorem. 

To estimate the tail of the sum over / in (p2D, we choose Iq so large that for / > /q, 
Ci{k) < 1/2 for all k in Nj. This is possible because the the compactness of a implies 
that its eigenvalues decrease to zero. Thus for each j we have 



l>lo 



l>lo^^^ 



Y^ / cnk)iPiu,(i),),,kdk 



< 



N. 



J2\iu,ipi)^'{ipi,(f)j)^'k\dk, (88) 



^ l>lo 



where the (pi are the normalized eigenfunctions of a. Here it may be necessary to 
reindex the sum. Each of the sequences {u, (pi)x' and [ipi, 4>j)x' is in P, and the inner 
product of two P sequences is in l^. We therefore find that the sum over / is bounded 
by llMllfcll^jlU, where || ■ ||a; denotes the norm in the space L|. Thus we have 



l>lo 



< 



N, 



\u\\k\\(l)j\\kkdk 



(89) 



which shows that the tail of the sum converges to zero as n goes to infinity. 

The same arguments, of course, apply to the denominator of (|8^). Thus we see 
that the leading order behavior of (^21) is given by the expression 



(t/„,$ 



B W 



Ei,jfi,jikj)J^^rJi{k)dk 
^i,j9ij{kj)Ur^{k)dk' 



(90) 
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where the sum in / is over those values for which A; attains the maximum M at kj. 
The terms ^^. rQ{k)dk, however, go to zero for large n. We must therefore consider 
their behavior in more detail. 

As we have seen, in the neighborhood oi k = kj, ro(fc) has an expansion of the 
form ro{k) = 1 — bj{k — kjY^ + . . ., where the positive integer pj is the order of the 
eigenfunction at kj. The lemma shows that the order pj controls the speed with which 
Jj^,rQ{k)dk goes to zero with n: the larger pj, the more slowly ^r^dk coverges to 
zero. We divide the numerator and denominator of (^) by /^ r^dk corresponding 
to the slowest decay. Finally, we take the limit of the resulting quotient as n goes to 
infinity. This shows that the quotient (|8^) converges to 

(TT ^ \ El,jPj{PluA)x'{kj) f f^ Piu{kj,x') 2 /,, 

(91) 
where the I3j are proportional to (3j = kjC{pj)/hj , and where the sums are over 
those indices j and / for which A/ has maximal order at kj. 
QED. 

Corollary 1 Assume that the sound speed c{x) is hounded and differs from cq only 
in a hounded suhset of the half-space x^ < —h < 0. Then for test functions \E'b, $b 
in X whose Fourier transforms with respect to time are supported in —B < k < B, 
Un as defined hy ^6^ ) converges to 

where the sums are over those indices j and I for which \i{kj) attains the maximum 
M and has maximal order. 

Proof. To apply Theorem 2, we need only check that a = J-'^^{VSV)*{VSV)J-' is 
indeed an analytic compact-operator-valued function of k. Analyticity was shown in 
section A.l; compactness was shown in section A. 3. The largest eigenvalue Aq cannot 
be constant: S is zero at fc = 0, which implies that all the Aj are zero there. Thus if 
Ao were constant it would be zero, and 5* itself would be zero. 



We note that both the numerator and denominator in (|92D can be zero, in which 



case (p2|) is not defined. However we do not study this case since the denominator is 



non-zero for a generic test function "^b- QED 
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C Appendix: Results from "Standard" Scattering 
Theory 



The solution of the equation 



{V^ + k^-V{x))ij{k,x) = 



(93) 



corresponding to an incident plane wave and a scattered field satisfying outgoing 
boundary conditions satisfies the Lippmann-Schwinger integral equation 



^jJ{k^ x,ri) = exp{ikri ■ x) — k I g{k^ x, y)V{y)^jj{k^ yiV)d- y 



(94) 



where g is the usual outgoing Green's function (^i]) . 

The initial difficulty with the Lippmann-Schwinger equation is that the incident 
field has infinite energy. This difficulty, however, can be circumvented by multiplying 
the whole equation by |l^(a;)|^/^ |T3]. This converts (p3) into 



where 



Cik, X, fj') = Coik, X, fi') + P / K{k, X - y)C{k, y, fi')d^y 



C{k,x,fj') 

Co{k,x,fi') 

Vi/2{y) 

K{x,y) 



\V{x)\'/'^{k,x,fi'), 
\V{x)\'/^e"'^-'', 
V{y)/\V{y)\'/' 
\V{x)\'/'g{k,x,y)Vr/2{y) 



(95) 



(96) 
(97) 
(98) 
(99) 
(100) 



When V has compact support, ( P5| ) is an integral equation on a bounded region. 
The kernel K is an Hilbert-Schmidt-valued function that is analytic in the entire 
complex fc-plane. By the analytic Fredholm theorem [^, the integral equation (p5|) 
is therefore uniquely solvable everywhere except at a discrete set of values of k, (the 
"exceptional points") and moreover the solution ^ is a meromorphic function of k 
with poles at these exceptional points. In addition, the arguments of [jl| |T^ show 
that the only possible real exceptional point is A; = 0. However, when A; = 0, we 
also have k = 0, and ( P3[ ) reduces to the equation C = |K|^/^. Thus the operator 
(/ — k^K)^^ is analytic in a neighborhood of the real A;-axis. 

This argument shows that for each k, ( = \V\^^'^ip is in L^. Then the quantity 
needed in section A. 3, namely |||V^|^''^ /5'V^'?/'||, can be rewritten as ||i^C|| = ||-^(-^ ~ 
-^)~^Co|| < ll-^IIIK-^ ~ -^)~^IIIICo||- Moreover, by explicit computation, we see that 
WCoik, -jf]')]] < ce"'^*'!'' I, where V is supported in the region y^ < —h < 0. 
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